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Theory for structure and bulk-modulus determination 
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Applied Physics, Chalmers University of Technology and Goteborg University, SE-41296 Gothenburg, Sweden 

(Dated: February 22, 2003) 

A new method for direct evaluation of both crystalline structure, bulk modulus Bq, and bulk- 
modulus pressure derivative B'q of solid materials with complex crystal structures is presented. 
The explicit and exact results presented here permit a multidimensional polynomial fit of the total 
energy as a function of all relevant structure parameters to simultaneously determine the equilibrium 
configuration and the elastic properties. The method allows for inclusion of general (internal) 
structure parameters, e.g., bond lengths and angles within the unit cell, on an equal footing with 
the unit-cell lattice parameters. The method is illustrated by the calculation of Bo and B'q for 
a few selected materials with multiple structure parameters for which data is obtained by using 
first-principles density functional theory. 

PACS numbers: 61.50.-f, 71.15.Nc, 62.20.Dc 
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I. INTRODUCTION 

Calculations of the bulk structure and the bulk elas- 
tic properties play an important role in the physics of 
condensed matter iSii^i^ Bulk calculations help us un- 
derstand, characterize, and predict mechanical properties 
of materials in our surroundings, under extreme condi- 
tions, as in geological formations and setting^ and for 
industrial applications^^ Crystalline materials come in 
many different structures and, in contrast to isotropic 
materials, the structure description of crystalline materi- 
als may in general need multiple lattice parameters and 
an atomic basis. In this paper we discuss how to de- 
termine the equilibrium structure of a (multiparameter) 
crystalline material while, at the same time, directly de- 
termining the bulk modulus and the bulk modulus pres- 
sure derivative. We argue and show that for theoretical 
structure calculations of multiparameter systems this is 
simpler and more exact than fitting to (semi-)empirical 
equations of state (EOS) such as, e.g, the Murnaghan or 
Birch EOS. In particular, with our direct method there 
is no need to first determine the hydrostatic path of the 
system. We further discuss how to include the atomic 
basis in this process in a natural way. 

In crystalline materials described by a single lattice pa- 
rameter {e.g, monatomic cubic phases) the lattice param- 
eter is a simple function of the unit-cell volume, and the 
equilibrium volume thus uniquely determines the equi- 
librium structure, i.e., the value of the lattice parameter. 
This is not the case when multiple lattice parameters 
characterize the system and a whole range of lattice- 
parameter values can form the same unit-cell volume. 
The equilibrium structure of the material must then be 
found by fitting and minimizing the free energy within 
the multidimensional space of lattice parametersiSi^ 
Relevant variables describing the atomic basis {e.g., bond 
lengths or binding angles) may be included among the 
parameters, and the full set of lattice parameters and in- 
ternal (atomic basis) parameters are collected into the 
vector X, scaled to dimensionless form. The volume of 
the unit cell V{x) depends in a simple way on the val- 



ues of the lattice parameters describing the unit cell, but 
not on the internal atomic configuration. Nevertheless, 
we here treat the external and internal parameters on an 
equal footing. 

From theory bulk calculations the total energy (per 
unit cell) E{x) is found for a number of structures x. 
The elastic response of typical hard crystalline materi- 
als corresponds to small deviations Sx = a; — aj*^"^ of 
the structural parameters from the equilibrium structure 
x^^\ The observation that the total energy forms a nat- 
ural potential (hyper-) surface in the parameter space of 
lattice and internal parameters x, combined with the ac- 
curacy of present-day bulk-calculation methods (such as 
density functional theory, embedded atom methods, or 
effective medium theory) , then makes it possible to fit the 
corresponding total-energy variation through the multi- 
dimensional fit 
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at controlled accuracy. Here k, M, and 7 denote zeroth, 
second, and third-rank tensors of fitting constants. An 
additional set of fitting constants are the a;'"' hidden in 
Sx — x — a;(°) The polynomial fit ^ gives a transpar- 
ent description of the materials-structure energy varia- 
tion and directly determines the equilibrium structure 

In this paper we exploit and use the structure cal- 
culation, i.e., the multidimensional polynomial fit of 
the total energy for an additional and direct de- 
termination of the zero-pressure bulk modulus Bq — 
~V{x'^°^){dp/dV)\^^^io} and its pressure derivative, 
Bq = —d{V{dp/dV)}/dp\^^^(o) at zero temperature. 

For a general set of structure parameters, x, we expand 
the volume around the equilibrium configuration a;'"' 
using the gradient g = VT^(a;)|^^^(o) and the Hessian 

H - H{V{x))\^^^,o, = \{d^V{x)/{dx^dx,)} ] 

of the volume. We note that derivatives of the volume 
with respect to the internal parameters vanish, by defini- 
tion. By providing a systematic treatment of the struc- 
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tural changes induced by the pressure p = —dE/dV we and the bulk- modulus pressure derivative 
extract from the minimum of the zero-temperature en- 
thalpy 

n{x,p)^ E{x)+pV{x) (2) 
both the bulk modulus 



B'o = V{x('y) 
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The algorithm outlined above can also be applied 
to the corresponding direct determinations of general 
harmonioifi and anharmonic elastic properties^ 

Our results both enhance the theory understanding of 
the crystalline mechanical properties and simplify the de- 
sired testing of theory calculations as they combine the 
formal determination of the crystalline structure [Eq. J^l] 
and of the elastic properties [Eqs. Q and Q]. For exam- 
ple, from Eqs. (QJ and (PI, we can directly identify which 
(internal) structure parameters softens the bulk modu- 
lus Q and we may, in turn, strengthen the materials by 
suitable chemical or structural modification. 

For cases like graphite, where the state-of-the-art DFT 
based on the generalized gradient approximation (GGA) 
fails to describe the weak physical interlayer binding, 
the polynomial fit ^ highlights the intrinsic theory 
challenges:^ the total energy has no minimum in the 
space of lattice parameters, the too soft dependence of 
the total energy on the graphite-layer separation is quite 
apparent in the fit. At the same time, for strongly bonded 
crystalline materials the fit Q produces usable esti- 
mates of the equilibrium crystalline structure and shows 
the strength of GGA-DFT. The equilibrium volume, the 
crystallographic parameters, and the bulk modulus, de- 
scribing the material's resistance to hydrostatic stress, 
provide simple experimental tests against which we can 
compare and calibrate our calculations. 

Besides the direct relevance of our results for the de- 
scription of complex materials our calculations of bulk 
structure and bulk modulus calculations are also of in- 
terest for development of pseudo-potential-based density- 
frmctional-theory (DFT) methods and for methods us- 
ing empirical parameters. There, a first and critical 
test of the pseudopotential or the empirical parameters 
is whether the calculations predict a correct materials 
structure, binding, and elastic properties for the relevant 
equilibrium configuration. Present DFT scripts^^ can au- 
tomate some pseudopotential testing for simple materials 
and symmetries, our formal results generalize such test- 



ing of theory accuracy to cases when multiple structural 
parameters determine the elastic properties. 

The outline of this paper is as follows. In Section II we 
discuss the traditional methods of determining the bulk 
modulus for single- and multiparameter systems. In Sec- 
tion III we derive our expressions for the bulk modulus 
and the bulk modulus derivative, Eqs. Q and for 
the simple one-parameter problem {e.g.^ mono-atomic fee 
or bcc structures), easily generalized to the n-parameter 
problem. In Section IV we proceed to illustrate and test 
the algorithm on a number of mono- and di-atomic mate- 
rials based on first-principle DFT calculations and com- 
parison to experiments. Comparisons of Bq and i?Q, to- 
gether with the test of the lattice and structure param- 
eters themselves, represent the typical test of materials- 
theory accuracy. Section V contains the conclusion. 



II. BACKGROUND 

A theory determination of the zero-temperature bulk 
modulus based on either traditional methodsitiiM or our 
formal result ||2Jl is straightforward when one single struc- 
tural parameter {e.g., the lattice parameter a) defines the 
crystalline state. This situation applies for monatomic 
crystals with simple cubic (sc), face-centered cubic (fee) 
and body-centered cubic (bcc) symmetries. Here, the 
unit-cell volume V{a) = qa?' uniquely determines the lat- 
tice parameter a through a dimensionless number q which 
depends on the crystal symmetry (g = 1, q = 1/4, and 
q — 1/2 for sc, fee, and bcc lattices, respectively). All 
which is required are theory calculations of total energies 
for a range of a values to determine both the equilibrium 
structure oq and the equilibrium volume V'^'^\ The total 
energy per unit cell, E{a) (as in Eq. (^l), can then be 
expressed as a function of the unit-cell volume, Eiy). 

The general approach is illustrated by the example in 
Figure n which shows the total energy as a function of 
the lattice parameter a for the zinc-blende phase of SiC 
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4.2 4.3 4.4 4.5 
lattice parameter a [A] 

FIG. 1: The total energy per unit cell (two atoms) as a func- 
tion of the lattice parameter a for the 3C polytype of SiC. The 
circles are the data points obtained from DFT calculations. 
Solid line: Fourth-order polynomial fit as used for the values 
in Table |I] Dashed line: Second-order polynomial fit to the 
15 central data points. 



(3C-SiC), as found from DFT calculations. A parabola 
fit to the 15 data points closest to the equilibrium value of 
a and a fourth-order polynomial fit to all data points are 
shown. Fits using the traditional Murnaghan equation 
of state/^ integrated to give 
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B'o-l 
(5) 



or the Birch equation of state, integrated to 
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^Bir(l^) = ~Eo + -BoV^°^ 
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(6) 



yield, to the eye, curves identical to the fourth-order poly- 
nomial fit and are not shown separately. In Murnaghan's 
and Birch's equations © and © the quantities Bq, B'q, 
and and in some cases also the cohesive energy Eq, 
are fitted. Other equations of state traditionally used are 
mentioned in Rcfs. Ill and 

The values of the equilibrium lattice parameter oq, and 
of Bq and Bq obtained from Eqs. Q and l@J and from 
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FIG. 2: Contour plot of the fourth-order polynomial fit to the 
Co hep total energy, including the hydrostatic path. The hep 
unit cell is given by the two independent lattice parameters 
a and c/a. The contour step is 0.025 eV per unit cell (two 
atoms) . 



the Murnaghan and Birch fits are included in Table 
For systems described by one lattice parameter Eqs. (O 
and lO give bulk moduli and bulk-modulus derivatives 
in close agreement with our present direct approach, 
Eqs. ® and Q. 

We would like to stress that the moduli Bq and Bq 
are formally defined as zero-pressure quantities, and in 
no way depend on finite-pressure behavior beyond the 
pressure gradient at p = 0. If we are able to sample 
our theory system in a sufficiently dense grid around the 
zero-pressure structure the values of Vq, Bq, and Bq in (O 
and Q are exact and can be related to the correspond- 
ing exact determination of the elastic constants. Fits to 
empirical EOS may yield results of Vb, Bq, and Bq that 
are in good agreement with experimental observations, 
but they do not necessarily constitute the exact quadratic 
response. 

For materials with multiple structure parameters, the 
procedure of the traditional approaches further becomes 
quite awkward as it must be supplemented by a separate 
discussion of how the experimental conditions define the 
relevant structural constraint at a given volume, the hy- 
drostatic path X = x{V). Moreover, cross-correlations 
on the rt-dimensional energy surface are ignored in tra- 
ditional fitting procedures. These procedures are basi- 
cally a one-dimensional fit in the n-dimensional space 
and they are thus more subject to numerical noise in the 
data points than our approaches based on the multidi- 
mensional least-squares polynomial fit 

A simple multi-parameter case illustrates this point. 
Figure 121 both describes the energy surface, Eq. (Q), fit- 
ted through DFT calculations for Co, and emphasizes 
the general advantages of our direct approach [Eqs. © 
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TABLE I: Bulk properties calculated from DFT data obtained directly from Eqs. @ and ||1J via a fourth-order polynomial 
fit ("Present approach"), available experimental values, and values from fits to Murnaghan's and Birch's equations Q and 
® along the hydrostatic path. For the internal parameter in 2H-SiC we find the following value: Si-C distance along the 
c-direction u(Si-C)= 0.3752c or bondlength it^nd = 1.9031 A. (Experiments findi^ u(Si-C)exp = 0.3760c, (4ond)exp = 1-8998 
A. ) For 4H-SiC we find li(Si-C)^ = 0.1880c, 'u(Si-C)2 = 0.1874c, ^(Si - Si) = 0.2500c, in good agreement with other 
theoretical resultsf^*— 
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and Ql] compared with the traditional bulk-modulus 
determinationsiiii^ Materials like Co, which has a non- 
ideal hexagonal close-packed (hep) structure, graphite 
with its layered structure, or the polytypes of SiG^ or 
alumina,^ have multiple lattice parameters (plus rele- 
vant internal degrees of freedom) which are, of course, no 
longer uniquely specified by the volume but depend on 
the general materials conditions. The hydrostatic path 
defines the system when subject to uniform pressure (as 
relevant for the bulk- modulus measurements). The hy- 
drostatic path is identified in Fig. [5] The traditional 
approachesi^*i^ proceed by implementing this complex 
constraint in the equation of state, Eq. (|3J) or ©, to fit 
the bulk modulus and its derivative. Instead, we present 
an explicit determination, Eqs. Q and l@J, based directly 
on the equation of state (Q) expressed as a function of the 
underlying crystalline structural parameters. 



III. DERIVATION 

Our direct bulk- modulus evaluations, Eqs. and 
are the results of using the pressure (instead of the vol- 
ume) in a formal identification of the hydrostatic path 
and then invoking a systematic expansion in small pres- 
sure for an explicit specification of Bq and B'q. To- 
day DFT and other materials-theory bulk-calculations 
are done with high accuracy, and we need only to vary 
the lattice parameters values slightly around the optimal 
structure to approximate the total-energy curve by an ac- 
curate polynomial fit, Eq. The minimum of the cor- 



responding (zero-temperature) enthalpy |(2J) can thus be 
used to directly specify the physically correct structural 
configuration at any given pressure p. The set of these 
optimal structure-parameter values, a^hydroCp), trace out 
the hydrostatic path which, when parameterized by p, is 
obtained by simply solving the equation 

Vnix,p) = VE{x) + pVV{x) = 0. (7) 

We obtain a formal expression for the general (pressure 
dependent) bulk modulus by taking the derivative along 
this hydrostatic path 

i?^-T/K...ob))(^^^^^^)"\ (8) 

and finally we extract the explicit results for the zero- 
pressure bulk-modulus Bq and for its pressure derivative, 
B'o- 

We illustrate the general derivation by focusing on a 
one-dimensional parameter space x, e.g., for the fee or 
bee one-atomic structure. Then the total energy can be 
fitted by the polynomial 

E{x) = k+-M{x-x^°'^f + -^-f{x-x^°^f + fix-x^°^) 

^ ^' (9) 

where f{x — x^^'') — 0{x — contains higher order 

terms. The coefficients k, M, 7, the coefficients of f{x — 
a;^°^), as well as the optimal value of the lattice parameter 
at zero pressure, x^^\ are the fitting parameters to be 
specified, for example, by a set of accurate underlying 
DFT calculations. 
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At small pressures, i.e., for lattice and internal pa- 
rameter values close to the zero-pressure optimal val- 
ues, the pV term in the enthalpy is small and can be 
regarded as a perturbation of the system. To proceed 
we introduce a small, non-dimensional and real param- 
eter A such that we can write the small pressure as 
p — Xp^^^ and the lattice-parameter variable as x — 
a;(o) + + A^cc^^) + 0{X^). The variables x^^), x'-^\ 
. . . are unknown, are functions of the pressure and must 
be found in the following. 

The bulk modulus expression requires calculation of 
the volume V{x) and its pressure derivative. We write 
the volume in a Taylor-expansion around the zero- 
pressure solution x^°^ as 



Vix) = F(x(")) + A.ga;(i) 



+AMx(2)g + 



H 



0{X^) (10) 



with g — dV/dx\^^^i_o) and H — d^V/dx^\^_^i^oy Here, 
the pressure dependence enters through the variables 
^ ^ _ _ The pressure derivative of the volume 
is thus 



dVjx) 
dp 



'dp(^) ^ 'dpW 



Hx 



(1) 



-0{X^) 

(11) 
by solving 



and we determine the variables x*-^-*, . 
the condition on the enthalpy given by ((TJ 

= A (a/x^i) + p(i),g) 

+A2 (^Afa;(2) -I- 1^ ' + p(i)i7x(i)^ + 0{X^) . 

(12) 

The identity p2|) must hold for every order and we thus 
obtain a formal pressure dependence of the lattice pa- 
rameter 



(13) 



,(2) _ 



(14) 



Finally, introducing these solutions into pi|l we find 
for A = the isothermal zero-pressure bulk modulus 



dp J x=xi») gM i<7 



(15) 



and taking the derivative of — (dV/dp) ^ with respect 
to p = Ap*^^) we find at A = 



^, ^ ^ 3gM-'HM-'g - ^M-^gM-^gM'^g 
° {gM-^gf 



- 1 
(16) 



in the case when one (lattice) parameter suffices to de- 
scribe the unit cell and its atom basis. 

The above derivation is straightforwardly generalized 
to materials systems in which n independent lattice and 
internal parameters determine the structure and the bulk 
moduli Eqs. and which is our main result. We 
stress that Bq and B'q are evaluated at zero pressure and 
thus the results are exact in spite of the perturbation. 
Bq and Bq depend directly on the second order, respec- 
tively on the second and third order, coefficients of the 
energy fit iM and 7). We observe that the coefficients 
of f{x — a;*-"-') do not enter the expression for the bulk 
modulus lj2Jl or the pressure derivative of the bulk modu- 
lus However, their presence may improve the fit 
and thereby affect also the coefficients A//, and 7, and 
thus _Bo and Bq. Internal parameters, which describe the 
positions of the atoms within the unit cell, naturally do 
not enter the expression of the volume, and thus not the 
volume derivatives g and H cither, but do affect Bq and 
B'q through A/-1. 

Higher pressure derivatives of the bulk modulus may 
be found by taking into account the higher orders of A in 
the Taylor expansions of x and the volume. The pressure 
derivatives will depend on successively higher orders in 
the polynomial fit. The derivation is straightforward if 
somewhat tedious. 



IV. EXAMPLES OF APPLICATIONS 

As an example of the use of the algorithm for deter- 
mining Bq and B'q we evaluate the structure and bulk 
modulus of a selection of one and two-species materials. 
We fit data obtained from DFT calculations, described in 
further detail below, to fourth-order polynomials of the 
form Q in n-dimensional space, where n = 1, 2, 3, or 5. 

The pseudopotentials used in DFT calculations may be 
optimized for various purposes, but should generally yield 
consistent and transferable accuracy and results. Here, 
we have used some of the pre-defined pseudopotentials^^ 
of the open-source DFT program DACAPO.^^ The values 
that we find for the lattice constants, for i?o, and for B'q, 
are collected in Table For reference, the experimental 
values are also included, as well as the bulk modulii from 
a Murnaghan and Birch fit along the hydrostatic path. 

We have calculated the structure and bulk modulus for 
three multiparameter systems, as well as for a number of 
related one-parameter systems: the two-parameter hep 
phase of Co, the three-parameter (one internal) wurtzite 
phase of SiC (2H), the five-parameter (three internal pa- 
rameters) hexagonal 4H-polytype of SiC, and the one- 
parameter bcc and fee phases of Co, the zinc-blende (3C) 
phase of SiC and the diamond phases of C and Si. 

For the DFT calculations we used the plane-wave 
DACAPO codei^ with GGA. For the calculations of the 2H 
and 4H-polytype of SiC we used 8x8x8 and 8x8x4 
fc-points, respectively, to describe the Brillouin zone. For 
all other calculations 10 x 10 x 10 fc-points were used. 
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A uniform energy cut-ofF of 400 eV, and a conservative 
choice of fast-Fourier transform (FFT) grid was used. For 
each evaluation of the optimal structure we calculated by 
DFT a number of data points for the lattice parameter(s) 
approximately within ±10% from the expected optimal 
value(s) of the lattice parameter(s). In the one-parameter 
systems we calculated 20-30 data points, for the Co hep- 
structure 120 data points, for the 2H-polytype of SiC 
140 data points, and for 4H-SiC we calculated 7500 data 
points. The Co calculations were spin-polarized, yielding 
realistic values of the spin polarization over the range of 
lattice parameter values considered here. 

The possibility of treating the external and internal 
parameters collectively is important. For example, the 
variation of internal bond length with pressure might be 
as important for the total energy as the change in (ex- 
ternal) lattice parameters. Often, relaxation of the inter- 
nal parameters is done with a steepest-descent (or simi- 
lar) search minimizing the HcUmann-Feynman forces to 
a certain cut-off at fixed lattice parameters. Although 
in practice the atomic relaxation will often be a conve- 
nient way of obtaining the optimal position of the atoms 
within the unit cell, the approach has two shortcomings: 
it introduces a random residual lattice strain, which in 
turn affects the total energy, and further, the Hellmann- 
Feynman forces have a non-trivial dependence on the 
pressure acting on the unit cell and therefore a constant 
cut-off on the force will not correspond to a constant ac- 
curacy of the total energy with varying pressure. Thus 
a better accuracy — and a consistent choice of accuracy 
— can be obtained by treating the lattice and internal 
parameters on an equal footing. This is here done for the 
2H and the 4H polytype of SiC. The results are shown in 
Tabled For the Murnaghan and Birch values of Bq and 



B'q we need to explicitly calculate the hydrostatic path 
[in (a, c/a, u) and (a, c/a, Ui, M2, M3) space] before obtain- 
ing the fit. In contrast, we stress that when using and 
Q there is no need to explicitly calculate the hydrostatic 
path. This is here done purely for illustrational purposes. 



V. CONCLUSION 

In summary, we have presented a new direct algorithm 
for a combined determination of structure and bulk mod- 
uli Bq and Bq. The lattice constants of a multiparameter 
system are best found in a least-squares polynomial fit, 
as previously noticed for the SiC hexagonal polytypesi^ 
We show (a) how to exploit this polynomial fit for a di- 
rect determination of the zero-pressure bulk modulus Bq 
and its pressure derivative Bq, avoiding the calculation of 
the hydrostatic path and the subsequent one-dimensional 
fit to this path. We further show (b) how to consis- 
tently include internal parameters, such as bond lengths 
or bonding angles, in the formalism along with the ex- 
ternal lattice parameters. In addition, we have evaluated 
these formal results in explicit cases within our approach, 
based on DFT calculations. 
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